• New Page
  • Convert a netcdf file to geojson

    The previously installed libraries must be:

    library("terra")
    library("geojsonio")

    Open several netcdf files and transform it to lat/lon projection and transform it to polygon. The files are located in the data folder:

    nc <- rast(c("./data/WRFDETAR_01H_20240101_12_000.nc",
                 "./data/WRFDETAR_01H_20240101_12_001.nc",
                 "./data/WRFDETAR_01H_20240101_12_002.nc",
                 "./data/WRFDETAR_01H_20240101_12_003.nc"))
    
    T2.positions <- which(names(nc) == "T2")
    
    T2.nc <- nc[[T2.positions]]
    
    T2.nc.repro <- project(T2.nc, "+proj=longlat +datum=WGS84")
    
    T2.polygons <- lapply(T2.nc, as.polygons)

    Finally, a shp file is written and the file_to_geojson function is used to create the geojson file for each time:

    for (i in 1:4){
      writeVector(T2.polygons[[i]], paste0("temps_", i, ".shp"), overwrite = TRUE)
    
      s <- file_to_geojson(paste0("temps_", i, ".shp"), method = "local",
                           output = paste0("temps_", i))
    }
    ## Success! File is at temps_1.geojson
    ## Success! File is at temps_2.geojson
    ## Success! File is at temps_3.geojson
    ## Success! File is at temps_4.geojson

    View geojson in Leaflet map

    The previously installed libraries must be:

    library("sf")
    library("leaflet")

    First the geojson to be graphed is opened:

    T2.geojson <- read_sf("temps_1.geojson") 

    Finally, the geojson variable is graphed, in this case temperature:

    pal <- colorNumeric("viridis", NULL)
    
    mapa <- leaflet(T2.geojson) %>%
      addTiles() %>%
      addPolygons(stroke = FALSE, smoothFactor = 0.8, fillOpacity = 0.7,
                  fillColor = ~pal(T2),
                  label = ~paste0(T2, ": ", formatC(T2, big.mark = ","))) %>%
      addLegend(pal = pal, values = ~T2, opacity = 1.0)
    
    mapa

    Ploting MCD12Q1 MODIS Product

    The data is available in data folder. First, open the data with terra package

    list.pathname <- list.files("./data", pattern = "hdf", full.names = TRUE)
    
    data <- lapply(list.pathname, FUN = rast)

    Now, creation of mosaic with each hdf.

    data.mosaic <- mosaic(data[[1]], data[[2]], data[[3]], data[[4]], data[[5]], data[[6]], data[[7]], data[[8]], data[[9]])
    ## 
    |---------|---------|---------|---------|
    =========================================
                                              

    Lastly, ploting of data with no change of projection, the numbers represent soil type

    plot(data.mosaic$LC_Type1)

    If the projection should be in longlat run:

    data.mosaic.repro <- project(data.mosaic$LC_Type1, "+proj=longlat +datum=WGS84", method = "bilinear")
    ## 
    |---------|---------|---------|---------|
    =========================================
                                              

    And plot again:

    plot(data.mosaic.repro)